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ABSTRACT 

A new theory of eccentric accretion discs is presented. Starting from the basic 
fluid-dynamical equations in three dimensions, I derive the fundamental set of one- 
dimensional equations that describe how the mass, angular momentum and eccentric- 
ity vector of a thin disc evolve as a result of internal stresses and external forcing. 
The analysis is asymptotically exact in the limit of a thin disc, and allows for slowly 
varying eccentricities of arbitrary magnitude. The theory is worked out in detail for 
a Maxwellian viscoelastic model of the turbulent stress in an accretion disc. This 
generalizes the conventional alpha viscosity model to account for the non-zero relax- 
ation time of the turbulence, and is physically motivated by a consideration of the 
nature of magnetohydrodynamic turbulence. It is confirmed that circular discs are 
typically viscously unstable to eccentric perturbations, as found by Lyubarskij, Post- 
nov & Prokhorov, if the conventional alpha viscosity model is adopted. However, the 
instability can usually be suppressed by introducing a sufficient relaxation time and/or 
bulk viscosity. It is then shown that an initially uniformly eccentric disc does not re- 
tain its eccentricity as had been suggested by previous analyses. The evolutionary 
equations should be useful in many applications, including understanding the origin 
of planetary eccentricities and testing theories of quasi-periodic oscillations in X-ray 
binaries. 

Key words: accretion, accretion discs - celestial mechanics ~ hydrodynamics - MHD 
- turbulence - waves. 



1 INTRODUCTION 
1.1 Background 

In a thin accretion disc the stresses due to collective effects 
are relatively weak and the motion of the gas is nearly bal- 
listic. Since the gravitational field is dominated by the cen- 
tral mass, the principal motion consists of Keplerian orbits. 
In the classical theory of accretion discs (e.g. Pringle 1981) 
these orbits are assumed to be circular and coplanar. Not 
only is this the simplest situation, it is also the expected 
outcome of the enhanced dissipation of energy that would 
result from initially misaligned or eccentric orbits. 

Nevertheless, there are strong observational and theo- 
retical grounds for believing that accretion discs in various 
situations are not flat but warped (e.g. Ogilvie 2000 and ref- 
erences therein). Such a situation can be caused by external 
forcing or can arise spontaneously through instabilities of an 
initially flat disc. Similarly, there are equally good reasons 
for investigating the possibility of eccentric discs. (The gen- 
eral case of a warped and eccentric disc deserves attention 
but is beyond the scope of the present paper.) 

The recent discovery of extrasolar planets orbiting 
main-sequence stars (Mayor & Queloz 1995; Marcy, Cochran 



& Mayor 2000) is a major development in the history of as- 
tronomy. The distribution of orbital elements of the plan- 
ets discovered to date poses some important challenges to 
theoretical models of planetary formation and dynamics. In 
particular, all of the planets with semi-major axes greater 
than 0.2 AU have significant eccentricities, several exceeding 
0.5. One of the leading contenders for driving eccentricity is 
the tidal interaction between the planet and the disc from 
which it forms (e.g. Lubow & Artymowicz 2000). 

A companion object on an eccentric orbit has a compli- 
cated tidal interaction with the disc. This is of considerable 
importance in connection with planetary rings and young 
binary systems, in addition to extrasolar planets. In the clas- 
sical theory (Goldreich & Tremaine 1980), tightly wrapped 
density waves are launched at an array of resonances and 
propagate some way through the disc before dissipating. 
The resulting torques lead to evolution of the orbit of the 
companion and of the surface density of the disc. In ad- 
dition, the companion may couple to global, low-frequency 
eccentric motions of the disc which are not described by the 
classical theory. The importance of such motions has been 
recognized by Tremaine (1998), and the methods required 
to understand them will be presented in this paper. 
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Further arguments support the notion of eccentric discs. 
Precessing eccentric discs are understood to exist in super- 
hump binaries as the result of a tidal instability (e.g. Lubow 
1991a and references therein). The general-relativistic apsi- 
dal precession of an eccentric inner disc is one of the leading 
explanations for kilohertz quasi-periodic oscillations in X- 
ray binaries (Stella, Vietri & Morsink 1999; Psaltis & Nor- 
man 2000; but see Markovic & Lamb 2000 for a critical as- 
sessment). A transitory eccentric disc may be produced after 
tidal disruption of a star close to a black hole at the cen- 
tre of a galaxy (Gurzadyan & Ozernoy 1979). Many-body 
systems with related dynamics include planetary rings, sev- 
eral of which have small but accurately measured eccentric- 
ities (e.g. Borderies, Goldreich & Tremaine 1983 and refer- 
ences therein), and the nucleus of the galaxy M31, for which 
Tremaine (1995) has proposed the model of an eccentric Ke- 
plerian disc. 

The existing theoretical work on eccentric discs con- 
sists of analytical studies and numerical simulations, based 
almost exclusively on a two-dimensional description of the 
disc. Numerous authors have treated small eccentric pertur- 
bations of an initially circular disc as a special case of wave 
modes. Kato (1983) showed that global, low- frequency ec- 
centric modes exist in an inviscid disc. A slow precession of 
the modes occurs owing to pressure forces. In a strictly invis- 
cid disc, eccentric modes exist that have non-trivial vertical 
structure (Okazaki & Kato 1985), although these are not 
expected to be important when viscosity is included. A sim- 
ilar, modal description was used by Hirose & Osaki (1993) 
to describe superhumps in SU UMa binaries. Ostriker, Shu 
& Adams (1992) considered the near-resonant excitation of 
eccentric density waves in a self-gravitating disc due to a 
companion on an eccentric orbit. Lee & Goodman (1999) 
used novel techniques to analyse non-linear eccentric den- 
sity waves in the tight- winding limit, for a self-gravitating 
disc. 

Much effort has been devoted to explaining the super- 
hump phenomenon in terms of a precessing, eccentric disc. 
Lubow (1991a) showed that a large-scale eccentric pertur- 
bation can grow through coupling with the tidal potential 
in a circular binary system. The coupling involves waves 
that are launched at eccentric Lindblad resonances, which 
are present only in tidally truncated discs in binaries of ex- 
treme mass ratio, as in SU UMa systems. Two-dimensional 
numerical simulations have been performed by Whitehurst 
(1988), Hirose & Osaki (1990), Whitehurst & King (1991), 
Lubow (1991b), Whitehurst (1994) and Murray (1996, 1998, 
2000). 

More general studies are of direct relevance to the 
present paper. Borderies et al. (1983) derived evolutionary 
equations for slightly eccentric planetary rings subject to 
quadrupolar and tidal forces, self-gravitation and viscosity. 
Syer & Clarke (1992, 1993) used the Gauss perturbation 
equations of celestial mechanics to argue that, contrary to 
the assumptions of classical accretion disc theory, an isolated 
viscous disc will retain its initial eccentricity indefinitely. 
Their numerical simulations supported this idea. Most re- 
cently, Lyubarskij, Postnov & Prokhorov (1994) derived evo- 
lutionary equations for a viscous accretion disc in which the 
orbits are ellipses sharing a common longitude of periastron. 
They confirmed the result of Syer & Clarke (1993) but only 
in the sense of a singular solution of the time-dependent 



equations, finding that an initially circular disc is viscously 
unstable to eccentric perturbations if the usual alpha viscos- 
ity model is adopted. 

The aims of the present paper are to unify these previ- 
ous treatments and to extend them in a number of important 
ways. A new set of evolutionary equations for an eccentric 
disc will be derived. The analysis of Lyubarskij et al. (1994) 
will be extended to allow for precession of the orbits, which 
can never be avoided and is the dominant feature in the 
pressure-driven modes of Kato (1983). The earlier modal 
description will be extended to allow for arbitrary eccen- 
tricities, and the effects of turbulent stresses and radiation. 
The analysis will also include important three-dimensional 
effects, which have been neglected in previous studies. In ad- 
dition, the conventional alpha viscosity model will be gen- 
eralized into a Maxwellian viscoelastic model, which is a 
physically motivated and more realistic description of the 
turbulent stress in an accretion disc. 



1.2 Continuum celestial mechanics 

As an introductory exercise, consider the problem of a test 
body orbiting in the gravitational field of a central object of 
mass M, but subject to a perturbing force per unit mass /. 
Its equation of motion is 

du GM 

where u = dr/dt is the velocity, with r = r Br being the 
position vector with respect to the central object. Since in- 
clined orbits are beyond the scope of this paper, assume that 
the orbit and the perturbing force lie in the a;i/-plane. 

The traditional method of determining the rate of 
change of the osculating orbital elements of the body in- 
volves evaluating the disturbing function and applying the 
Gauss perturbation equations of celestial mechanics (e.g. 
Brouwer & Clemence 1961). A more compact derivation uses 
the eccentricity vector e, which lies in the plane of the or- 
bit and may be defined by (Eggleton, Kiseleva & Hut 1998; 
Lynden-BeU 2000) 



GM 

u = -J— X (Br + e), 

or, equivalently, 
h 



GM 
where 

he^ = r X u 



■ e, X U — Br 



(2) 



(3) 



(4) 



is the specific angular momentum. The position and veloc- 
ity are then instantaneously equal to those of an elliptical 
Keplerian orbit of semi-latus rectum 



^=GM=' + '-^ 



(5) 



and eccentricity e — |ej, with e pointing towardsperiastron. 

Substitution into the equation of motion {}^ gives the 
relations 



dh 

■^e.-rxf, 

de dh , N , j: 
'^dt = dt ^^'■ + ^)~^^^^^' 



(6) 
(7) 
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which determine the evolution of the orbital elements and 
are equivalent to the Gauss perturbation equations. 

The problem to be addressed in this paper is how to 
derive a related set of equations for a fluid disc in which 
e = e(A, t) varies not only with time but also from one orbit 
to the next. The evolutionary equations will be partial differ- 
ential equations (PDEs) depending on one spatial variable, 
analogous to those obtained for warped discs (Pringle 1992; 
Ogilvie 1999, 2000), and equally suitable for semi-analytical 
or numerical implementation. 

The perturbing force in a fluid disc includes a contri- 
bution from the internal stress in addition to any external 
forcing. The instantaneous elliptical Keplerian motion as- 
sociated with the function e{\,t) determines (part of) the 
velocity gradient tensor in the disc. An appropriate theory, 
satisfying the basic principles of continuum mechanics, is 
then required to relate the stress tensor to the velocity gra- 
dient tensor. 



1.3 Modelling the turbulent stress 

The simplest assumption, made in almost all theoretical 
work on accretion discs, is that the turbulent stress may 
be treated as an isotropic effective viscosity in the sense of 
the Navier-Stokes equation.^ The stress tensor may then be 
written as 

T = ^l\^u + {Vu)'^]+{^l^-l^i){v■u)l, (8) 

where ^ is the effective (dynamic) shear viscosity and /^b 
the effective bulk viscosity (which is often omitted) . For the 
usual case of a steady, flat and circular Keplerian disc, this 
assumption suffices to generate the stress component Th^ 
that is required for accretion to proceed. Together with the 
conventional parametrization of the viscosity coefficient, this 
constitutes the alpha viscosity model (Shakura & Sunyaev 
1973), which was very successful in connecting theory and 
observations during a period when the theory was incom- 
plete. 

In recent years it has become widely accepted that 
the turbulent stress in most, if not all, accretion discs is 
magnetohydrodynamic (MHD) in origin, resulting from the 
non-linear development of the magnetorotational instability 
(Balbus & Hawley 1998). The turbulent stress tensor is 



B'B' 

4tv 



II 1 1 I , I , II 
pu u — p (uu + u u + u u 



(9) 



where the prime denotes a turbulent fluctuation from the 
mean value of a quantity, and the magnetic pressure fluc- 
tuation has been considered to be combined with the gas 
pressure fluctuation. Numerical simulations of the turbu- 
lence offer the possibility of measuring the stress and test- 
ing the validity of the alpha viscosity model in more general 
circumstances. In fact very few studies have been directed 
towards this aim (Abramowicz, Brandenburg & Lasota 1996; 
Torkelsson et al. 2000). 

The most obvious deficiency of the alpha viscosity 
model is that it requires the stress to respond instanta- 
neously to a change in the rate of strain. In an eccentric 

* Note the crucial distinction between isotropy of the stress and 
isotropy of the effective viscosity. 



disc the strain field experienced by a fluid element changes 
during the course of each orbit. There is a strong suspicion 
that the assumption of an instantaneous response may be 
partially responsible for the eccentric instability described 
by Lyubarskij et al. (1994). In reality the turbulence is ex- 
pected to have a non-zero relaxation time comparable (on 
dimensional grounds) to the orbital period. 

It is therefore proposed that the stress tensor satisfies 
the equation 

T + rVT = /i [Vu + {Vuf] + (^b - |/i)(V-u)l, (10) 

where r is the relaxation time and V a linear operator de- 
fined by 

VT = dtT + U-VT - (Vu)'^ ■ T - T-Vu + 2(V-it)T. (11) 

This model (although with V-u = for an incompressible 
fluid) is used in rheology as a simple model of a viscoelastic 
medium, known as the upper- convected Maxwell fluid. The 
relaxation time was introduced by Clerk Maxwell (1866) in 
his study of the viscosity of gases. The Maxwellian viscoelas- 
tic fluid allows for a continuum of behaviour from an elas- 
tic solid (on time-scales short compared to r) to a viscous 
fluid (on time-scales long compared to r). The operator "D is 
one of a family of convective derivative operators, acting on 
second-rank tensors, that satisfy material frame indifference, 
which is a fundamental principle of continuum mechanics. In 
the context of rheology, equation ( |lo| ) can be justified either 
by abstract principles of continuum mechanics, or from a ki- 
netic theory of polymer molecules that satisfy Hooke's law 
(Bird, Armstrong & Hassager 1987; Bird, Curtiss & Arm- 
strong 1987). 

There are good reasons to suppose that a viscoelastic 
model may offer a reasonable description of MHD turbu- 
lence in accretion discs. Magnetic field lines have tension 
and support Alfven waves similar to waves on an elastic 
string. There is a close analogy between the tangled mag- 
netic field lines in a turbulent medium and the tangled poly- 
mer molecules in a viscoelastic fiuid. Indeed, the equation 
satisfied by the dyadic tensor BB in ideal MHD is precisely^ 

V{BB) = 0, (12) 

and {B' B' / in) is arguably the most important contribution 
to the total turbulent stress tensor (^). 

This accounts for the 'elastic' aspect of MHD turbu- 
lence. The 'viscous' part of equation (|l^) can then be related 
informally to the (less well understood) dissipative and non- 
linear aspects of the turbulence. Viscoelastic models have 
been proposed for hydrodynamic turbulence (Crow 1968) 
but the case for MHD turbulence is arguably stronger. 

This hypothesis is also related to the 'causal viscosity' 
model used in the context of accretion disc boundary layers 
to ensure that information is not propagated faster than the 
sound speed (Kley & Papaloizou 1997). In that context, only 
a single component of T is treated and the full expression 
for V is not used. Here, however, it is essential to treat the 
stress as a tensor and to ensure that the equations satisfy 
material frame indifference. 

Note that one easily recovers the conventional alpha 
viscosity model by taking the limit r ^ 0. 

t I am indebted to Prof. M. R. E. Proctor for pointing this out. 
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It is of interest to work out the implications of this 
model for a steady, flat and circular Keplerian disc. The 
horizontal stress components are predicted to be 



Trr = 0, 



where Q is the angular velocity and 



(13) 



(14) 



the Weissenberg number, which is expected to be of or- 
der unity. A comparison with the results of MHD simula- 
tions by Hawley, Gammie & Balbus (1995) suggests that 
the Maxwellian viscoelastic model performs rather well in 
predicting the turbulent stress tensor (^). The results of Ta- 
ble 2 in that paper imply that 



{Trr , Tr^ , Tr2 , Ttf,^ ,T42,Tzz) 

(X (0.105, -1, 0.018, 1.520, -0.050, 0.076), 

whereas the viscoelastic model predicts 



(15) 



{Trr,Tr4„Tr,,T4,4„T4,,,T,,) (x (0, -1, 0, 3We, 0, 0). (16) 

Given that the numerical quantities are averaged over a 
limited time interval and involve some statistical error, 
the agreement is good and the results are consistent with 
We ~ 0.5. It would be valuable to conduct further simu- 
lations designed to measure the relaxation time in a more 
direct way. 



1.4 Plan of the paper 

The remainder of this paper is organized as follows. In Sec- 
tion 2 I present a highly simplified analysis of an eccentric 
disc, with the purpose of explaining the principles of the 
method without the technical detail of the full model. In 
Section 3 I set out the basic equations and coordinate sys- 
tems used in the rest of the paper. Section 4 contains the 
bulk of the analysis. I present an asymptotic development 
of the equations appropriate for a thin disc, reducing the 
problem to simpler units that are solved in turn to extract 
the required one-dimensional evolutionary equations. In Sec- 
tion 5 a linearized theory for small eccentricity is described 
and compared with the oversimple model of Section 2. An 
illustrative, time-dependent solution of the fully non- linear 
equations is given in Section 6. Finally, in Section 7, the re- 
sults are summarized and a comparison with existing work 
is made. Recommendations for applications of this work are 
also given. 



Suppose that the fluid is strictly two-dimensional and 
obeys the equation of mass conservation. 



{dt + u-V)p = —p\/-u, 

and the equation of motion, 

p{dt + u-V)u = -pV$ - Vp -f V-T, 



(17) 



(18) 



where p is the density, u the velocity, p the pressure and T 
the turbulent stress tensor, which will be assumed to sat- 
isfy equation (p^. The gravitational potential of the central 
mass M is $ = —GM/R, referred to plane polar coordinates 
(i?, 0). The self-gravitation of the fluid is neglected. 

Much use has been made of the analogy between the 
dynamics of a two-dimensional fluid and that of a thin disc. 
However, the analogy is not a formal one except under spe- 
cial circumstances (and never for an eccentric disc, as will 
be seen). In fact, a two-dimensional model corresponds to 
an inflnite cylinder rather than a thin disc, although the 
gravitational potential cannot be interpreted in this sense. 
However, p, p and in the two-dimensional theory might 
be regarded as approximately equivalent to the correspond- 
ing vertically integrated quantities in the three-dimensional 
theory. To emphasize this, in this section the density will 
be written as E (surface density) instead of p. To avoid a 
consideration of the energy equation at this stage, suppose 
that the fluid is barotropic, so that p — p(S), p — /i(E) and 
Mb = /^b(S) are given functions. 

Let {u, v) be the polar components of u. Then one has 



dt + udR + ^a^) ^ = [dR{Ru) + d4,v] , 



(19) 



{dt + uOr + -^9^,^ ^~~R 
+ ^dR{RTRR) + i^^Tij^ 



'^-^-dnp 



R 



dt + uOr + —d^y ) " + 
+ ^dR{R''TR^) + ^d^T^^ 



(20) 



(21) 



with 



Trr + t^ (dt + udR -f- ^d4^ Trr + ;|('^ + d4,v)TRR 
-{d4,u)TR^'^ = 2pdRU 



2 



+{p^ - |m)-^ [dR{Ru) + d^yv] , 



(22) 



2 AN OVERSIMPLE MODEL 
2.1 Analysis 

In the following sections a detailed model of a thin, eccen- 
tric accretion disc will be developed, consisting of a fully 
non-linear asymptotic analysis of the fluid-dynamical equa- 
tions in three dimensions, including dissipation of energy 
and radiative transport. Before setting out this theory, I 
present a grossly simplifled model (two-dimensional, linear 
and barotropic) that is not recommended for further use 
but nevertheless illustrates some of the important principles 
without many of the technical details of the full calculation. 



■[{dt + udR + ^a^) Tr^ - RdR (^) 

i [dR{Ru) + d^yv] Tr^ - ^{d4,u)T^^^ 



Trr 



RdR ( ^ + ^d^u 



(23) 



T^^ + r I (dt + udR + ^d4,^ T4,4, - 27?^^ (^^) Tr^ 
+2idRu)T^^} = + ia^^) 



-K(^ib - [dR{Ru) + d^yV] 



(24) 
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Let the basic state be an axisymmetric, circular disc 
that evolves on the slow, viscous time-scale. Adopt a system 
of units such that the radius of the disc and the character- 
istic orbital time-scale are 0(1). The slow evolution is cap- 
tured by a slow time coordinate T = e^t, where the small 
parameter e is a characteristic value of the angular semi- 
thickness H/R of the disc. For the basic state, introduce the 
expansions 

u = e\2{R,T) + 0{e^), (25) 
V = Rn{R) + €\2iR,T) + 0(e*), (26) 
S = Eo(i?,r) + O(e'), (27) 
p = 6^ [po(i?,T) + 0(e')] , (28) 
M = [/io(i?,r)-KO(e')] , (29) 

Mb = [/ibo(ii,T) + 0(e')] , (30) 
T = [To(7?,r) + O(e')] , (31) 

where 



1/2 



(32) 



is the Keplerian angular velocity. Substitution into equations 
(^-(|l|) leads to 



(9t -I- u2i9i{)Eo = -—dR{Ru2), 



2Eofi«2 = -Orpo — 



R ' 



1 1 2 



with 

Trro — 0, 
Tfl^O ~ — 7 



9 



(33) 
(34) 
(35) 

(36) 
(37) 

(38) 



and this leads to the well known evolutionary equation for 
the surface density (e.g. Pringle 1981), 



(39) 



Now consider slowly varying, one-armed linear pertur- 
bations of the basic state, such that the Eulerian perturba- 
tion of u, say, is 

Re[u'{R,T)e-''^] . (40) 
For the perturbations, introduce the expansions 

= u'o{R,T)+e^u'2{R,T) + 0{e*), (41) 
= v'o{R,T) + t\'2iR,T) + 0{e*), (42) 
= K{R,T) + 0{e^), (43) 
= [p'o{R,T) + Oie')] , (44) 

m' = [/io(i?,r) + 0(6')] , (45) 

Mb = [f^UR,T) + 0{e^)] , (46) 

T' = e"" [To{R,T) + 0{e^)] . (47) 

In this linear analysis, the overall amplitude of the perturba- 
tions is of course arbitrary, and the above scaling is chosen 
for convenience. 

Equations (|^) and (|2l]) at 0(1) yield 



Eo{~mu'o - 2nv'o) = 0, (48) 

Eo(-i!^^;o + ^nu'o) = 0, (49) 
with the general solution 

u'o = iRQE, (50) 

v'o = ^RnE, (51) 

where E{R, T) is a complex function to be determined. It 
is easily verified that this solution corresponds to a small 
eccentric perturbation, with complex eccentricity E — e e'" 
corresponding to eccentricity e(R,T) and longitude of peri- 
astron u>{R,T). Equivalently, E = -I- ifiy, where {ex,ey) 
are the Cartesian components of the eccentricity vector. 

Equation ( p^ ) at 0(1) then yields the density pertur- 
bation, 



= RdR(EoE). 



(52) 



The pressure and viscosity perturbations follow from the 
barotropic relations, e.g. 



(53) 



Equations ( |2c| ) and ( ^l| ) at O(e^) contain u'2 and v'2, 
but these may be eliminated by taking an appropriate linear 
combination. This results in an evolutionary equation for E, 
which may be written in the form 

2T.oR^dTE = -2EoU2dR{RnE) - QdR{R'£oU2E) 



-2iR^'^mR{R^'^Y.oV2E) + ±-dR{R:'p'a) 

R-' 



-^dR{RT'R, 



2R 



-'"•''^dRiR^^Tr 



R4,Q) 



R ■ 



The stress perturbations satisfy the equations 

Tflflo + We{-iTitRo - 2Tr^oE) 

= 2iiiodR{RnE) + i(/ibo ~ liio)RndRE, 



T'r^o + We 



Trro — VE'r^q -f iRdR{TR^oE) — T^^qE 



^\ (Is) + ^^^o^RiR'ilE), 



r;^o + We [sn^o - ir^^o - R'^^dR{R-'^^E)TR^Q 
+iREdRT^4,o + 2iR'^^dR{R-^'^E)T^^^^ 
= ifio^E + i(/ibo - ifio)R^dRE. 



(54) 



(55) 



(56) 



(57) 



Thus one has obtained evolutionary equations for the 
surface density (equation ^9|) and the complex eccentricity 
variable (equation . In this linear theory the equations are 
decoupled because the small eccentricity makes a negligible 
change to the stresses that cause accretion. 



2.2 Evolution of eccentricity in a purely viscous 
fluid 

Equation (|5^) may be written in the form of a linear evolu- 
tionary equation for the eccentricity. 



2T,oR^dTE = AR^OrE + BROrE + CE. 



(58) 
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The general form of the complex coefficients (.4, B, C) is com- 
plicated, but in the limit of a purely viscous disc, We 0, 
one obtains 



f ^^^\ , 2,,^ , ^ / dp 

Q VdE 



(59) 



+ [(/^^ - + -^d~{R^^%) , (60) 

C = -2RdRii + '-§dRp, (61) 

where the subscript zeros have been omitted. 

First, suppose that an initial condition is specified in 
which E = constant, so that the eccentricity is uniform and 
the ellipses are aligned. The initial rate of change of E is 
determined by the coefficient C. In a steady disc far from 
the inner radius, is nearly independent of 7? and the real 
part of C nearly vanishes. Therefore there is no initial vis- 
cous evolution of eccentricity, as has been argued by Syer & 
Clarke (1992, 1993) and Lyubarskij et al. (1994). However, 
the imaginary part of C does not vanish (except in the un- 
likely case that the pressure is independent of R) and causes 
a differential precession of the orbits. Viscosity will then act 
on the twisted configuration. It follows that a steady, uni- 
formly eccentric disc is not a solution of the model. This 
difference comes about because the calculations of Syer & 
Clarke and Lyubarskij et al. neglected the differential pre- 
cession caused by slightly non-Keplerian rotation resulting 
from the radial pressure gradient. 

Now consider a limit in which the eccentricity varies 
on a length-scale £ <^ R. (One requires H, however, in 
order that the ordering scheme adopted above remain valid.) 
Then equation (|54[) becomes approximately 



drE ^ — 



1 

2E 



3/1 



m) 



iE 



/dp 
n VdE 



diE. 



(62) 



This is a linear diffusion equation with a complex diffusion 
coefficient. The imaginary part of the coefficient, which is 
due to pressure, results in wave-like propagation of the ec- 
centricity and to the modes studied by Kato (1983). The real 
part determines the viscous diffusion of eccentricity. If the 
real part of the diffusion coefficient is negative, instability 
occurs. Within this oversimple model, therefore, an initially 
circular disc is unstable to developing eccentricity on small 
scales if 



dln/i ^ , 1 f Mb _ 2 
dlnE 3 I /i 3 



(63) 



i.e. if the vertically integrated viscosity increases sufficiently 
rapidly with increasing surface density. This criterion agrees 
with the findings of Lyubarskij et al. (1994) if the bulk term 
(the term in brackets) is neglected. However, this analysis 
shows that a sufficiently large bulk viscosity stabilizes the 
disc. 

Interestingly, equation (^) is identical to the criterion 
for viscous overstability of short-wavelength axisymmetric 
modes in a two-dimensional disc (Kato 1978; see Willerding 
1998 for the effect of bulk viscosity). This should not be 



surprising because the local properties of waves in a thin disc 
are almost independent of the azimuthal wavenumber when 
the radial wavelength is much shorter than the azimuthal 
wavelength. A related result was also found by Goldreich 
& Tremaine (1978), who examined the damping effect of 
shear and bulk viscosity on density waves in planetary rings. 
They considered the case in which jj, and /Xb were constant, 
but commented that viscosity might produce growth rather 
than decay if /i and /ib were to depend on E. Note that the 
condition (|6^) is not related to the criterion for thermal- 
viscous instability, dln/i/dlnE < 0. 

In the alpha models (Shakura & Sunyaev 1973) for gas- 
pressure dominated discs, din /j./ din is equal to 5/3 for a 
Thomson opacity law and 10/7 for a Kramers opacity law. 
Stability therefore requires /ib/M to exceed 8/3 or 41/21, re- 
spectively. Although bulk viscosity is not often discussed in 
the context of alpha models, there is no reason to suppose 
that accretion disc turbulence has a vanishing effective bulk 
viscosity. In addition, it is known that radiation damping, 
which has been neglected here, can mimic a bulk viscosity 
(Papaloizou & Pringle 1977). If the instability does occur, 
however, it will produce short-wavelength variations that in- 
validate the theory developed here. The non-linear outcome 
would have to be followed using numerical simulations that 
allowed for the viscosity to be a function of the surface den- 
sity. The non-linear development of the axisymmetric ver- 
sion of the instability has been computed by, e.g., Papaloizou 
& Stanley (1986) and Kley, Papaloizou & Lin (1993). 

Starting from equation ( |54| ) it is possible to derive a 
conservation law for eccentricity in the form^ 
1 



dT{T.R^fl\E 



2R 

-|-viscous terms. 



^%R''{iE*dRE - lEdRE*) 

QZj 



(64) 



This shows that, in the absence of viscosity, the eccentric- 
ity is conserved in a certain sense. The conserved quan- 
tity is related to the 'angular momentum deficit' which is 
conserved in the secular theory of celestial mechanics (e.g. 
Laskar 1997). The eccentricity flux due to pressure effects 
is similar to the particle flux in Schrodinger's equation. The 
viscous terms are not necessarily negative deflnite, indicat- 
ing the possibility of instability as described above. 



2.3 Effect of a non-zero relaxation time 

The general form of the coefficient A for a Maxwellian vis- 
coelastic fluid is 



A = a 



-iWe)-^{3M-3E(^) +(Mb-iM) 
(t^) [(l + iWe)M + (Mb-|M)]} 



3i We 
iWe 

Q. VdEy 



(65) 



Instability on small scales occurs if Re(^) < 0, and one finds 
(1 + WefReiA) = 3(1 + me - We*))j, 



"3( 



l + We2)Ef4^U(l 



VdE/ 



(1 + 7We') {pb - (66) 



t I am indebted to Dr S. H. Lubow for suggesting this calculation. 
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For either Thomson or Kramers opacity, the instability per- 
sists for all values of We if /^b/M = 0. However, for reasonable 
values of the bulk viscosity and relaxation time, the insta- 
bility can be quenched. For example, the case ^ih/ ^J^ ~ 2/3 
corresponds to the vanishing of the 'bulk term' in equa- 
tion (10|). For this value, the instability is quenched for 
3-1/2 < -y^g < 2^/2 (i.e. 0.577 < We < 1.414) in the case of 
Thomson opacity, and for 0.423 < We < 1.547 in the case of 
Kramers opacity. It is tantalizing that the value tentatively 
deduced from simulations of MHD turbulence. We « 0.5, 
may or may not be sufficient to stabilize the disc, and the 
outcome may differ according to the opacity law. 

A related result is that the viscous overstability of ax- 
isymmetric waves can be quenched by introducing a more 
sophisticated turbulence model with a non-zero relaxation 
time (Kato 1994). However, a direct comparison with Kato's 
work is not possible. His second-order closure model is aimed 
at modelling the turbulent Reynolds stress in an incompress- 
ible fluid that undergoes a notional hydrodynamic instabil- 
ity, whereas the Maxwellian viscoelastic model used in the 
present paper is aimed at modelling (predominantly) the 
turbulent Maxwell stress in a compressible fluid that under- 
goes the magnetorotational instability. 



3 THE FULL MODEL 
3.1 Basic equations 

In this section the basic equations governing a fluid disc 
in three dimensions are expressed in a general, time- 
independent coordinate system. Following the usual nota- 
tion of tensor calculus, gab denotes the metric tensor and 
Va the covariant derivative. 

The equation of mass conservation is 



{dt + U'Va)p = -pVaU", 



(67) 



where p is the density and u the velocity. The equation of 
motion is 

p{dt + u'VtK = -PV"* - V> -I- VoT"" + pr, (68) 

where $ is the gravitational potential, p the pressure, T the 
turbulent stress tensor and / the external force per unit 
mass. 

Under the assumption that the turbulent stress behaves 
similarly to a Maxwellian viscoelastic fluid, the stress tensor 
satisfies equation (|ic|), i.e. 

T"'' + T [{dt + uVc)T'''' ~ T'"'VcU - T'"'V^u -f 2(Vcm")T'''' 

= pi^^u'' + V'u") + (Mb - I m)(V.^*= )<;"', (69) 

where r is the relaxation time, p the effective shear viscosity 
and /ib the effective bulk viscosity. The energy equation is 



{dt 4- u°'Va)p = 



7 



7— ly \7— 1 



(70) 



where 7 is the adiabatic exponent and i^rad the radiative 
energy flux, given in the Rosseland approximation for an 
optically thick medium by 

3 

(71) 



where a is the Stefan-Boltzmann constant, T the temper- 
ature and K the opacity. The equation of state of an ideal 
gas, 



P = 



kpT 
Pmirm ' 



(72) 



is adopted, where k is Boltzmann's constant, /im the mean 
molecular weight and mn the mass of the hydrogen atom. 
The opacity is assumed to be of the power-law form 



(73) 



where Ck is a constant. This includes the important cases 
of Thomson scattering opacity (x = 3/ = 0) and Kramers 
opacity {x = 1, y = —7/2). 

The gravitational potential is $ = —GM/r, where r 
is the distance from the central object. The self-gravitation 
of the disc is neglected. The effective viscosity coefficients 
are assumed to be given by an alpha parametrization. The 
precise form of the alpha prescription relevant to an eccentric 
disc is to some extent debatable. It will be convenient to 
adopt the form 



p — ap 



Ph = QbP 



(74) 



where a and Qb are the dimensionless shear and bulk vis- 
cosity parameters, and A is the semi-latus rectum of the 
elliptical orbits (introduced below). In the limit of a circular 
disc, this prescription reduces to the usual one, /i = ap/Q,, 
etc., where Q is the orbital angular velocity. Finally, the 
Weissenberg number is defined by 



We = r( — j 



{75) 



where r is the relaxation time. Again, this reduces to We — 
Q,T in the limit of a circular disc. 



3.2 Orbital coordinates 

Consider a thin, eccentric accretion disc in which the domi- 
nant motion consists of coplanar, elliptical Keplerian orbits. 
Let {R, 4>, z) be cylindrical polar coordinates such that the 
mid-plane of the disc corresponds to 2 = 0. Then the shape 
of the orbits is given by 



i? = A(l-l-ecose')' 



(76) 



Sup 



where A is the semi-latus rectum, e the eccentricity and = 
<j) — uj the azimuth relative to periastron (known as the 'true 
anomaly' in celestial mechanics). The semi-latus rectum is 
related to the specific angular momentum h — {GMXf^^ 
and is conveniently adopted as a label for the diff'erent orbits. 
In general the eccentricity and longitude of periastron may 
vary continuously from one orbit to the next and also with 
time, so that e = e{\,t) and u = u!{X,t). This is illustrated 
in Fig. 1. Let e and e' then denote dte and dxe, etc. It 
will also be convenient to define the complex eccentricity 
E = ee'"'. 

The analysis of an eccentric disc is greatly facilitated 
by the introduction of orbital coordinates (A, (p) in place of 
{R,4>). However, some preliminary analysis must be carried 
out because these coordinates are both non-orthogonal and 
time-dependent. 
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Figure 1. Example of instantaneous orbits in an eccentric disc. 



3.2.1 Spatial derivatives 

To evaluate the spatial derivatives in the governing equa- 
tions, one must obtain the metric coefficients and connec- 
tion components. The two-dimensional line element is given 

by 

ds^ = AR^ + d(^^ = {Rx dA + R^ d0)^ + R^ d(^!>^ (77) 

where the subscripts on R stand for partial derivatives of the 
function R{\,(j>,t). These can be evaluated from equation 
([tgI), but the more general notation will be retained. Thus 
the metric coefficients are 



<7aa — R\, g\4, — R\R, 



■4> 5 



^R^ + Rl. (78) 



As expected, the square root of the metric determinant is 
equal to the Jacobian of the coordinate system; 



1/2 ^ J ^ d{x,y) _ d{x,y) d{R, 



d{X,<j)) d{R,cl>)d{X 
The inverse metric coefficients are 



= RRx. 



9 



R^Rl ' 



R6 



(79) 



(80) 



The components of the metric connection, given by 

rjc = ^g'"^{dbgcd + dcQdb - ddgtc), (81) 

are found to be 



Rx 



R 



\<t> 



\<t> 



R(t> 



rtA = o 



{R^ + 2i?0 — RR^ 
RMx 
^ _ Rx 



(82) 

(83) 
(84) 



p0 _ 

A0 - ^ ' 4"t>- R ■ 

The coordinate system is trivially extended to three di- 
mensions by incorporating the vertical coordinate z. With 
the exception of gzz ~ g"" ~ 1, all metric coefficients and 
connection components involving z vanish. The divergences 



of a vector and of a symmetric tensor in three dimensions 
are then 



Va ^ a , T^a b 

aU = da It -I- i baU 



bT = dbT + TcfcT + TcfcT . 



(85) 
(86) 



The vector divergence can be written explicitly as 

jdx{Ju^) + jd^iJu"^) + dzu\ (87) 

while the A-, 0- and z-components of the tensor divergence 
are 



1 / ^AA^ , R'' a ( JR>^rpX4> 



dci, 



i?2 



^ (i?' + 2Rl ~ RR4,4,)T'^'>' + dzT^^ , 



(88) 
(89) 
(90) 



RR 

-l^dxiJR^T^'^) + --^d4,(JR^Tl"f)+dzT'^ 

jdxiJT^n + jd^iJT^n + 
respectively. 

For a physically meaningful solution, the eccentricity is 
restricted in magnitude. In order that the orbits be closed, 
one obviously requires \E\ < 1. A further restriction is that 
the Jacobian determinant should remain positive, which im- 
plies \E — XE'\ < 1. If this condition were violated, neigh- 
bouring orbits would intersect each other. 



3.2.2 Time-derivatives 



The equations of Section 3.1 are valid only in a time- 



independent ('inertial') coordinate system. To adapt them 
to orbital coordinates = {X,4),z), observe that the par- 
tial derivatives at fixed inertial and orbital coordinates are 
related by 



(9t)inGrtial — (9t)orbital + q'^da, 

where 

q" = (9t)inortiaig". 



(91) 



(92) 



Thus, for a scalar quantity such as p in equation (^^ or p 
in equation ()70|), one replaces 



dtp dtp + q'dap. 



(93) 



When differentiating the vector u in equation (|68|) one must 
also take into account the time-dependence of the orbital 
basis vectors. One should then replace 

dtu^ ^ dtu^ + qdbu' — u'dbcf . (94) 

This is most easily derived by noting that = u-VA, etc., 
and obtaining the equation for the time-derivative of this 
quantity. More generally, one is replacing dt with dt + 
where £ denotes the Lie derivative. Thus in equation 
one should replace 



dtT" 



(95) 



Hereafter dt will represent the derivative at fixed orbital 
coordinates. For orbital coordinates only A is non-zero, so 
if da = A^A. The identity 



dtJ + dx{JX) 



(96) 
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will be used below. This is established by noting that 

(9t)inertial-R = 0, SO that 

dtR + \Rx = 0, (97) 
and so 

dtJ + 9a(JA) = dt{RRx) - dx{RdtR) = 0, (98) 
as required. 

4 ASYMPTOTIC DEVELOPMENT 

4.1 Scaled variables and expansions 

The key to analysing this problem, as with many dynamical 
problems in accretion discs, is to recognize that there is a 
separation of scales. The disc is thin, in the sense that the 
semi-thickness H of the disc satisfies H/R <^ 1. Addition- 
ally, the evolution of surface density and eccentricity occurs 
on a time-scale much longer than the orbital time-scale of 
the fluid. These ideas lead naturally to the introduction of 
scaled variables and an asymptotic development. 

Let the small parameter e be a characteristic value of 
the angular semi-thickness H/R of the disc. Adopt a system 
of units such that the radius of the disc and the character- 
istic orbital time-scale are 0(1). Then define the stretched 
vertical coordinate ( — z/e, which is 0(1) inside the disc. 
The slow evolution is captured by a slow time coordinate 
T — e^t and one may replace 

e(A,t) e(A,r), Lj{\,t) lo{X,T). (99) 
For the fiuid variables, introduce the expansions 



= eV(A,<^,C,r) + 0(e*), 



(100) 



T 



u 

= n{\,^,T) + e\ti\<P,(,T) + Oie^), (101) 

= euU\,c^,(:,T) + e\liX,<t>X,T) + 0{e''), (102) 

p = e= [po(A,<^,C,r) + e2p2(A,0,C,T)-f 0(e*)] ,(103) 

= ^+'[T,^\X,<t>X,T) + 0{e^)], 



r 



s+2 



e-'-[T,^^iX,<P,(:,T)+0{e^)] 

[ri^^(A,0,c,r) + o(e^)] 

T"^" = [rf"(A,0,C,T) + O(e'')] 



Mb = 

r = 



[ro"(A,<^,C,T) + 0(62)] 
[^*o(A,0,C,T) + O(62)] , 
^o(A,<^,C,T) + 0(e')] , 

[ro(A,<^,C,r) + o(62)] , 
;j'o(A,<^,C,r) + o(6')] , 



s+3 



(104) 
(105) 
(106) 
(107) 
(108) 
(109) 
(110) 
(111) 
(112) 
(113) 
(114) 



while the horizontal components of i^rad are ©(e^"*"*). Here 
s is a positive parameter, which drops out of the analysis, 
although formally one requires s = (4 — 2y) /(2 -I- a;) in order 
to balance powers of e in the opacity law. Note that the 
dominant motion is an orbital motion with angular velocity 
SI, independent of 

The potential is expanded in a Taylor series. 



where 



$0 = - 



GM 



$2 



GM 



R ' ^ i?3 ■ 
The external force per unit mass is also expanded as 

= e^[fo\X,^,T) + 0{e^)], 
f = e^[f^{XA,T) + 0{^)], 

r 



(116) 

(117) 
(118) 
(119) 



This is the correct form if / represents the tidal force due to 
a companion object in a coplanar orbit. The scaling will be 
seen to be appropriate, and indicates that the external force 
is small compared to the gravitational force of the central 
object. 

When these expansions are substituted into the equa- 
tions of Section 2, various equations are obtained at different 
orders in e. The required equations comprise three sets, rep- 
resenting three different physical problems, and these will 
be considered in turn. 



4.2 Orbital motion 

The horizontal components of the equation of motion 
at leading order [0(e°)] are 



(120) 



po{nd^^ + r^^fi') = -po(<7^^aA$o + p'^^^^^o). (121) 

Since <l>o is a function of R only, these simplify to 



{R^ + 2R^ — RR^^)Q^ — R ,^ , 
d^R^^l) = 0, 



(122) 
(123) 



and are satisfied by the angular velocity of an elliptical Ke- 
plerian orbit. 



(124) 



where e{X,T) and u!{X,T) are arbitrary. It is useful to note 
from the above that 



R^ + 2i?0 — RR^tf) — , 



(125) 



which simplifies the expression for one of the connection 
components. 



' XRx' 



The orbital period is^ 



1/2 



(126) 



(127) 



$ = <l>o(A, T) + ie\'$2(A, 0, T) + 0(e*), 



(115) 



S Throughout, integrals with respect to <f> are carried out from 
to 27r, and integrals with respect to C, over the full vertical extent 
of the disc. 
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4.3 Slow velocities and time-evolution 

The equation of mass conservation at leading order 

[0(.= )] is 



{0.84, +uldc)po ^ ~po 



— d^{JQ) + d(ul 



(128) 



Introduce the surface density at leading order [©(e""*"^)], 
E(A,0,T) = / podC. (129) 



Then the vertically integrated version of equation ( 12J ) may 
be written in the form 



d^iJtn) = 0, 



(130) 



which determines the variation of surface density around the 
orbit due to eccentricity. It will be convenient to introduce 
a pseudo-circular surface density, 



E(A,T) = — / JE, 



(131) 



which has the property that the mass contained between 
orbits Ai and A2 is 



27r / EAdA. 
/a. 



(132) 



It will also be useful to introduce the second vertical moment 
of the density, 

i(A,0,T) = J poC'dC, (133) 

and the pseudo-circular average, 

I(A,T) = ^ / 

The equation of mass conservation ( |67[ ) at 0(e''^^) is 
[Ot + {u2 + \)d\ + U2d4, + uld(] po + (0,84, + uldc)p2 



-po 



-p2 



(135) 



The horizontal components of the equation of motion 
at 0(e"+^) are 

- g oxpo - g ^O4>po 



2Rx dR 

+;^a.(ji..To-) + -^a,(:^T- 
-i^rr + acT'^ + Po/o, 



Po [dr + (u2 + A)9a + utdity] ^ 

+P0 [i^id4, + uldc)ut + 2rt^nu^ + 2r^^n^i^] 



(136) 



g ^oxpo ~ 



+ —dxiJRX'") 



+ -±^d4JR^T,^^) + dcTt + Po/o^ 



(137) 



It will be sufficient to wo rk with the vertically integrated 
versions of equations (135), (136) and (137), which may be 



written in the form 

[Ot + \dx)t+jdx j Jp^u\ dC 

^-\^^4. J J{pout + p2n)dC = 0, (138) 

nd^iv^ -\) + 2r^x4.^v^ + 2r^^nv'^ = a^, (139) 

[Ot + {v"- + \)dx + v'^d4,] ^ + ^d^v"^ 

+2rt^nv' + 2^^^n^;^ 

where t;(A, (j), T) is a mean horizontal velocity defined by 
tv"" = / pouadC, (141) 



(140) 



and j4(A,(j!),T) a mean horizontal acceleration given by 



EA^ = -J^^ + -^dxURxT'') 



+ 



2Rx dR JRx 



JRl 



i?2 



A-T^^ + E/o\ 



(142) 



tA^ = -^dx{JR^T^^) + -^d^iJR^T-^^) + tf^, (143) 
where 



(144) 



is the vertically integrated stress tensor including the pres- 
sure term. 

The objective at this stage is to extract evolutionary 
equations for E and E without attempting to obtain a com- 
plete solution for v'^ and v"^. First, equation (138) is inte- 
grated over to eliminate the mass fiux around the orbit 
and obtain 



dr / Jtd(t) + dx / Jt{v^ + X)d(t> = 0, 



(145) 



where the identity (^6|) has been used. This expresses the 
conservation of mass in one dimension and involves the rel- 
ative velocity [v^ + A). Let v{X,T) be a mean accretion ve- 
locity defined by 



V jt.d(P= JE(u^ + A) 



(146) 



Then equation (145) may be written in the form 



(147) 



E+-9;,(A,;E) =0, 

formally identical to the equation of mass conservation for 
a circular disc. 



Next, equations (139) and ( |140| ) may be rewritten in 
the form 



A " R^ ' 



(148) 



(149) 



nd^{Rv'')+v''— = R-'B 
qA 

where h = {GMXy^^ is the specific angular momentum and 

(150) 



= A^ + Qd4.X, 



i?2 dA 
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Equation (L49) expresses the conservation of angular mo- 



mentum and may be integrated over the orbit to obtain 



^ / JE(^;^ + A)d<^ = 



JE7?^yl"* d(j) 



(151) 



which involv es ex actly the same average relative velocity as 
in equation (145). Indeed, these two equations clearly relate 
the evolution of mass to the total torque, as expected on 
general grounds. This may be simplified to 



dh 



(152) 



dx p J n 

Using the expression for A'*, this becomes 

J Jtd(t> = dx j JR^T^''' d(t) + j jtR^f^dcj}. (153) 

It remains to determine the evolution of the eccentric- 
ity. Here one uses the fa ct th at th e op erator acting on the 
unknown v in equations (148) and (149) is singular, and the 
equations can therefore be solved (in principle) only when 
B satisfies a certain solvability condition. Eliminate v'^ be- 
tween equations (148) and (149) to obtain 

Cv'>' 



n, 

where the linear operator £ is defined by 



dA A 



(154) 



(155) 



(156) 



and the right-hand side is 

TZ = nd,iRlB^)-'^^B\ 

Now £ is a singular operator because it possesses a null 
eigenvector 



e 
R^ 



(157) 



such that Lw'^ = 0. This can be verifi ed by direct substitu- 
tion, with the help of equations (123) and (125). This null 
solution corresponds simply to a small redefinition of the 
eccentricity, E E + SE . The corresponding solvability 
condition for equation (154) is 



J'SR^w'^Tldtj) = 0. 



(158) 



Multiplying this by 2i\ /GM and taking into account the 
contributions to B from A, this can be brought into the 
form 



E I JEd 

2iA 

GM 



2^0 



JE 



Rx 



nd,{RlA^)-f^^A' 



where the relation 

-RaA e cos 6 + eu> sin 6 

R I + e cos 6 

has been used, and also the integrals 



(1 + ecos6')"''d6l = (2 -f e^)7r(l - e^) 



2N-5/2 



d(^, (159) 



(160) 



(161) 



(l + ecos6')"^cos0d0 = -3e7r(l - e' 



2N-5/2 



(162) 



(H-ecos6')~^cos^e'd6' = (1 2e'')7r(l - e^)"'^''^(163) 



which are easily established by the residue theorem. Thus 
the required equation for the evolution of E has been ob- 
tained, but must now be manipulated into a more usable 
form. 

With the help of the relation 



-RT^^y-RT 

one obtains 
E / Jtd<j) 



= i(e'^ + E- \E') 



(164) 



GM 



+ - 



JT.n{e"^ + E- \E')R''A'^ d<^. 



GMX 

This can be further simplified to 



(165) 



h[E + vE' - ^ \ = ^ I {2R^A'^ - iXRxA^) e'"^ ^-(166) 



X ) P 



Taking into account the various contributions to A, this 
becomes 

h[E + vE' I JEd0=-| I jm^e''^d(l} 



-iXdx j JRxT^^e"^d(f> + 2dx j JR^T^"* e"^ d(l> 
-i j JR^'T^'^ie"^ + E-XE')d<t> 



i / JR^'T'^'^e"*', 



+ / JE(2i?V(?-iAi?A/o)e''^' 



(167) 



Thus a complete set of one-dimensional equations has been 
obtained that determine the evolution of the mass, angular 
momentum and eccentricity vector due to internal stresses 
and external forcing. It remains to evaluate the equations 
explicitly for the case of a Maxwellian viscoelastic model 
with an alpha viscosity. 



4.4 Vertical structure and vertical motion 

The shear tensor of the orbital motion will be required be- 
low. This is defined in general by 



b =2(Vu-fVM), 
and has the expansion 
S<^^ = St + 0{e). 

The horizontal components at leading order are 
x\ {R^Rxtp -f RRxtpR^ -f RxR^ — RRxRtpRtpip)^ 



5o' = 



R3RI 



[R^ + Rl)dxO. + {RxR4,4, - R4,Rx4,)0. 



2R^Rl 



(168) 
(169) 

(170) 
(171) 
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■ {Rd^n + Rxn), 



(172) 



where equation (123) has been used. The orbital contribu- 



tion to the divergence Va^t" is 



(173) 



where gi{9;X,T) is a dimensionless quantity, equal to zero 
for a circular disc (see Appendix A). Note, however, that 
there is also a vertical shear component Sq^ — d^ul of the 
same order, which remains to be evaluated. This contributes 
to both the divergence and the dissipation rate at leading 
order. 

Noting that 

1/2 



2n~3/2 / ^'^ 
" \GM 

is the orbital period, one may write 
S = E;?2, 



(174) 



(175) 

where g2{9;X,T) is a second dimensionless quantity, equal 
to unity for a circular disc, and satisfying 

(1 + ecos6)^9£) In (72 = — (?i. (176) 

Further useful quantities are defined according to 



93, 



V A3 y 



94, 



and 



bo — 



(177) 
(178) 

(179) 



V A3 y 

so that (33, 514, s^'*', As^"^, A^s*^**) are dimensionless. These 
quantities are written out in Appendix A. 

The equation of mass conservation (^) at leading order 
[0{en] is 

{nd^, + uldc)po = -po(A + d^ul). (180) 
The energy equation (^^ at leading order [0(6"^^)] is 



7-1 



[0,84, + uld^)po 



po(A + d^ul) 



(181) 



The vertical component of the equation of motion (|68|) at 
leading order [0(e^"''^)] is 



(182) 



The required components of the stress equation ( |69| ) at lead- 
ing order [0(e°"'"^)] are 



T^^ + T [{nd^. + uldc)T^^ + 2(A + dcnDT,^"-] 
= 2/ioS'o'^ + {pho - f Mo)(A + d^uDg^^, 

+2(A-hc»c"i)ro'"^] 



(183) 



(184) 



-f2(A + 9c«i)ro^1 

= 2moS'o^^ + (Mbo - |mo)(A + 9c^i)5^^ 

To"" + r [(f7c»^ + uldc)T,r + 2Aro"] 

= 2fiodi;ul + (/ibo - |mo)(A + d^ul). 

The constitutive relations at leading order are 



Fo 



Po 



kpoTo 



^0 = apo 



/GMV 
V A3 ) 



-1/2 



MbO = «bPO 



/GMV 
V A3 ) 



-1/2 



(185) 
(186) 

(187) 
(188) 
(189) 



The solution of equations ( 18C )-(189) is found by a non- 
linear separation of variables, a similar method to that used 
for warped discs (Ogilvie 2000). As an intermediate step, 
one proposes that the solution should satisfy the generalized 
vertical equilibrium relations 



dpo 



-h po [— ) C, 



dFo , 9a fGM\^^^ 



(190) 
(191) 



in addition to equations (187)-(18£). Here fi{9;X,T) and 
f2(d; A, T) are dimensionless functions to be determined, and 
which are equal to unity for a circular disc. Physically, fi dif- 
fers from unity in an eccentric disc because of the enhanced 
dissipation of energy and because of compressive heating 
and cooling (equation 181). Similarly, /2 reflects changes in 
the usual hydrostatic vertical equilibrium resulting from the 
vertical velocity, the vertical viscous stress and the variation 
of th e ve rtical oscillation frequency around the orbit (equa- 
tion 182). Following Ogilvie (2000), one solves these equa- 
tions by reducing them to a standard dimensionless form. 
One identifies a natural physical unit for the thickness of 
the disc. 



Uh 



9ay/(«+--2'') (2+.)/(6+.-2«) (GM\ 



(t) 



-(5-2y)/2(6+x- 



PmrUH 



V A3 y 

-(i-y)/{a+x-2y) ^ •,-l/(6 + x-2y) 



and for the other variables according to 



Uo 



Ut 



Uf 



(GM\ //immn 
V A3 J V k 



(192) 

(193) 
(194) 

(195) 
(196) 



The solution of the generalized vertical equilibrium equa- 
tions is then 



^ ^ j.l/(6 + x-2y) ^-{3-y)/(li+x-2y) 
^^i2+x)/(e+x-2y)^^^ 



(197) 
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''Up, 



^g2(2-B)/(6+x-2i;), 



PO =P*(C*)/l 



l/(6 + a:-2a) .{3+a;-a) /{6+a;-2a) 
J2 



2(4+x-y)/(6+x-2y) 
^92 



Up, 



2(2+x)/{6+x-2y) 



•2/{6+x~2y) rx/(6+x~2y) 
■12 



X52 



Ut, 



(8+x-2y)/(6+x-2y) rx/(6 + x~2y) 
J2 



(198) 
(199) 
(200) 
(201) 
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(211) 



Fo = F.(C.)/{ 

y^giU)+3x-2y)/(&+x-2y)jj^ 

where the starred variables satisfy the dimensionless ODEs 
dp* 

d-F, 

dC* 

dT* 

p. 



= P* 



subject to tlie boundary conditions 

F.(o) = p.(Cs.) = r.(Cs,) = 

and the normalization of surface density, 



P.dc* = i. 



(202) 

(203) 

(204) 
(205) 

(206) 
(207) 



a. 



Here C,b* is the dimensionless height of the upper surface of 
the disc. Also required is the dimensionless second vertical 
moment of the density, 

= J P.C^dC.. (208) 

The solution of these dimensionless ODEs depends only on 
the opacity law and is easily obtained numerically (Ogilvie 
2000)0 For Thomson opacity one finds Cs. « 2.383 and 
X* = 0.6777. For Kramers opacity one finds (s* ~ 2.543 and 
J* = 0.7094. 

One further proposes that the vertical velocity is of the 

form 



c, 



(209) 



where faiO; \,T) is a third dimensionless function, equal to 
zero for a circular disc, and that the stress components at 
leading order are of the form 

TS' = t''''{e;\,T)po, (210) 

so that {t^^,\t^'l',\^t'l"^,f) are dimensionless coefficients. 
When these tentative solutions are substituted into 



equations (180)-(186) one obtains a number of dimension- 
less ODEs in 9, which must be satisfied if the solution is to 



^ It is assumed here that the disc is highly optically thick so 
that 'zero boundary conditions' are adequate. The corrections for 
finite optical depth are described by Ogilvie (2000). 



be valid. From equation ( |180[ ) one obtains 

(1 + e cos of [-de In /i + (3 - y)dB In /a] 
= -{6 + x-2y)h-{2 + x)gi. 



Similarly, from equation (181) one obtains 
1 



(1 + ecos f 



7-1 



de In /2 = 



7 + 1 
7-1 



h -91 



+t^\xx + 2t^^SA0 + t^^s^^ + t"/3 - ^a/i. (212) 



From equation ( 182 ) one obtains 

(1 + ecosOfdeh = -.fl - (1 + ecosef + /2(1 - f"").(213) 



Finally, from equations (183)-(186) one obtains 
t^^ + We [(1 + ecosef{d0t^^ + t^^de In /a) + (3/3 + <7i)t^^] 
= 2as^^ + (ab - |a)(gi + fi)g^\ (214) 



At^"* + We [(1 + ecos6')^(ae(At^"*) + At^'^ae In/a) 



+ (3/3+3i)At'^-33t^^-.g4At^^l 



= 2a As^"^ + (ab - |q)(3i + h)\9 

\^t1"^ + We [(1 + e cos ef{de{\^t't"l') + A^t^'^ae In /a) 
+ (3/3 + gi)\h^^ - 2g3\t^^ - 234 A^t^^] 
= 2a\^s'^^ + (ab - |a)(ffi + h)\^g'^\ 



(215) 



(216) 



+ We [(1 + e cos e)\dee^ + e^de In /a) + (/3 + gi)^^] 

= 2a/3 + (ab-fa)(5i+/3). (217) 

These ODEs should be solved for the functions 
(/i, /2, fs.,t^^, At^**, A^t''"*, t") subject to periodic boundary 
conditions fi{2n; X,T) — fi{0; X,T), etc. Note that, in the 
case We — 0, the equations for t"'' become algebraic, reduc- 
ing the problem to third order. Note also that, in the limit 
of a circular disc (e — 0, gi — 0, g2 — 1, gs — — |, g4 — 0, 
g^"- = 1, Ag^"* = 0, AV^ = 1, sxx = 0, \-^sx^ = -|, 
A~^S0^ = 0), the solution is fi = 1, /2 = 1, = 0, t^^ = 0, 
Xt^"^ = -fa, AV"^ = fWea, t"" = 0. 



4.5 Evaluation of the stress integrals 

The evolutionary equations may be written in the form 



t + jdxiXvE) = 0, 

^"dA = A^n^^^— J + p/ "^^"V' 



.{e + vE'-^4)=o.(z.i^]+Z2i'''' 



A 



A2 J 



A3 



where Qi is a dimensionless real coefficient defined by 



(218) 
(219) 

(220) 
(221) 



and Zi and Z2 are dimensionless complex coefficients defined 
by 
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i / Jii^T^'^e'*' d<^ = Z2 



GM 



(223) 



To evaluate the stress integrals, first note the relation 

(224) 



which follows from equation (19C) after an integration by 
parts. Then 



/T-ab /,ab ab\ r (GM\~ 

T ^{t -g )h[—]^- 

Define the dimensionless function /4(S; A,r) by 

Jl = /4AX, 

so that 

where 

/5 = (l + ecos0)-Vi' 



2 f2/{6+x-2y) j.-2(3-y)/(6 + x-2y) 



^2(2+x)/((i+x-2y) 

Then one finds 

(3i = ^ / /2/4(l + ecose)-'A(t^*'-5^'^)d6l, 



Z, = -^e'" / /2/4 [2(1 + ecose)-^A(t^^ - <?^^) 



(225) 
(226) 
(227) 

(228) 
(229) 

(230) 



+ 2^e'"y" |-|/4(l + ecose)'* + /2,f4(l + ecos6l)-" 

X [A(f^^ - 3^"*) - iA''(f'*'^ - 5^"^)] } e'" d6». (231) 
Finally, a relation is required between E and T. One 



finds 



2/(6 + a:-2H) / 'j^" \ 



I = Q2Cxa . . 

\ A-^ / 

^ j.(10 + 3a:-2i/)/(6 + a:-2i/) 



(232) 



where 



_ /9\2/(6+— 2^) /^i^nmH\-2(4-!/)/(6+==- 



2y) 



/ 16a \ -2/(6+— 2!/) 



V3ay 



(233) 



is a constant and Q2 is a dimensionless coefficient defined 
by 



Q2 = ^{l~e^f'^ / /sde. 



(234) 



and equal to unity for a circular disc. The dimensionless co- 
efficients can all be evaluated numerically from the solution 
of the dimensionless ODEs in Section 4.4. For given values 
of the model parameters (a, Qb, We, 7, a;, j/), the quantities 
(Qi, Q2, .^i/e"^, .Z2/e'") depend only on the dimensionless 
dynamical variables (e, Ae', Xeuj'). 



4.6 Relation to the Gauss perturbation equations 

It may now be shown how the evolutionary equations for 
a continuous disc relate to those derived in Section 1.2 for 
a test body. The relation between the ordinary cylindrical 
polar components (./^,/^) of the perturbing force and the 

is 



contravariant orbital components (/"^ , 

Thus equation (^) for the test body becomes 

dh dh 2 ,0 
d^^^dA^^^ ' 



(235) 
(236) 

(237) 



where v = dX/dt. This is equivalent to equation (152) for the 



disc, except that the disc equation involves a time-average 
over the orbit, and the perturbing force includes contribu- 
tions from internal stresses. 

Similarly, equation (^) for the test body may be (after 
some algebra) brought into the form 



(238) 



which is consistent, in the same sense, with equation (166) 



4.7 Relation to standard accretion disc theory 



I n th e absence of an external torque, equations (218) and 
( [219| ) may be combined into a single evolutionary equation 
for the surface density. 



-d.iX 



X 



.1/2. 



. /GMy/2 



(239) 



In the limit of a circular disc, A reduces to R and one finds 
that Qi — 3a/2. One then obtains the standard diffusion 
equation for the surface density of a circular accretion disc 
(e.g. Pringle 1981), 



(240) 



where 1/ — aQ{X/T,) is the vertically averaged kinematic 
viscosity. In an eccentric disc, therefore, the surface density 
diffuses viscously in much the same way as in a circular 
disc, except that it is best considered as diffusing in the 
space of semi-latus rectum A (which is equivalent to specific 
angular momentum). Non-linear couplings arise, however, 
because the coefficient Qi depends on the local eccentricity 
distribution, specifically on e, Ae' and Xeu'. The elliptical 
distortion of the orbits modifies the stresses that leads to 



accretion. The eccentricity obeys its own equation (220) 
which has the character of a complex, non-linear advection- 
diffusion equation. 
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5 LINEAR THEORY 

In general the dimensionless ODEs must be solved numer- 
ically and the Q- and ^-coefhcients determined from the 
solution. However, for small eccentricities the equations can 
be solved by expanding in powers of e. In doing so it is as- 
sumed that XE' is 0(e), i.e. that the length-scale on which 
e or w varies is comparable to the radius of the disc. 
One then finds 



Qa = 1 + 0(6^), 

Zi = ciE + C2XE' + 0{e^), 

Z2 = C3E + C4\E' + 0{e^). 



(241) 

(242) 
(243) 
(244) 



The complex coefficients Ci depend only on the dimension- 
less parameters (a, ab, We, 7, i, y), and can be obtained by 
a straightforward algebraic calculation. Such a calculation 
is not presented here because the resulting expressions are 
too complicated to write down, and it is no more difficult 
to evaluate the dimensionless coefficients numerically in the 
non-linear case. The only simple, general statement that can 
be made is that the Ci become purely imaginary in the in- 
viscid limit a, Q!b — > 0. 

If the higher-order terms are neglected, and in the 
absence of external forces, one obtains a complex linear 
diffusion-type equation for the eccentricity, together with 
the standard evolutionary equation for the surface density, 
which is decoupled from the eccentricity. This is similar to 
the situation in Section 2 except that three-dimensional and 
radiative effects have been included. Of greatest interest 
is the coefficient C2, because the small-scale instability dis- 
cussed in Section 2 exists if Re(c2) < 0. 

For the purposes of illustration, suppose that a = 0.1 
and 7 = 5/3. In a purely viscous disc with We = 0, it is then 
found that stability requires a bulk viscosity ab > 0.350 in 
the case of Thomson opacity, or Qb > 0.176 in the case of 
Kramers opacity. It is by no means clear that such values are 
realistic. However, the instability is easily quenched when a 
non-zero relaxation time is introduced. Even when ab = 0, 
the instability is absent for 0.398 < We < 2.294 (Thom- 
son) or 0.333 < We < 2.467 (Kramers). This contrasts with 
the oversimple model of Section 2, which generally predicted 
narrower stability intervals and suggested that the instabil- 
ity persisted for all We when /ib = 0. These results empha- 
size the importance of including three-dimensional efi'ects 
and radiative damping when treating eccentric discs. 

One of the features of the linear analysis is that shows 
that vertical motion is always present in an eccentric disc. 
To understand this, consider a simplified example in which 
the gas is isothermal so that the pressure is proportional 
to the density. Suppose further that the eccentricity is in- 
stantaneously uniform. In the absence of vertical motion the 
velocity divergence is zero and the density and pressure must 
then be constant along any orbit, for any value of z. How- 
ever, this is incompatible with hydrostatic vertical equilib- 
rium because the vertical gravitational acceleration varies 
around any eccentric orbit. Of course, this efi'ect is missing 
in two-dimensional treatments of eccentric discs. 



6 AN ILLUSTRATIVE NON-LINEAR 
CALCULATION 

The practical nature of the non-linear evolutionary equa- 
tions is now demonstrated by a simple illustrative calcula- 
tion. In any detailed application, careful consideration must 
be given to the formulation of appropriate boundary con- 
ditions and to any sources of mass and/or eccentricity to 
be added to the equatio ns. A t a free boundary of the disc, 
where E — > 0, equation (220) generally has a singular point 
and one must select the regular solution. If the disc is termi- 
nated by an external agency, however, a different boundary 
condition on the eccentricity may apply. These issues should 
be addressed in future work, within the context of specific 
applications. For the present purposes, it is convenient to 
study a simple test problem with idealized boundary condi- 
tions. 

The parameters adopted are a = 0.1, ftb = 0, We = 0.5, 

7 — 5/3, X = 1 and y = —7/2, appropriate to a disc 
in which Kramers opacity is dominant. (The correspond- 
ing linear-theory coefficients are ci = 0.075 -I- 0.972i, C2 = 
0.030+0.6311, C3 = -0.2284-0. 955i and C4 = -0. 198-0.2221. 
Therefore the short-wavelength instability is absent.) For 
this opacity law there exists a self-similar solution with 
E oc A"^'''' and E = 0, representing a steady, circular accre- 
tion disc having an arbitrarily small inner radius. According 
to the theories of Syer & Clarke (1992, 1993) and Lyubarskij 
et al. (1994), if such a disc is made uniformly eccentric, with 
all the orbits aligned, it should remain so. In order to con- 
trast the results of the present theory with these earlier cal- 
culations, the non-linear evolutionary equations (218)-(22C) 
are solved starting from an initial condition with surface 
density S oc A~^'^^ and uniform eccentricity E = 0.3. The 
equations are solved in a finite domain, 1 < A < 100, with 
boundary conditions 91nE/i91nA = —3/4 and dxE = at 
the inner and outer edges. These illustrative boundary con- 
ditions are chosen to be relatively neutral while being com- 
patible with the notional solution proposed by Syer & Clarke 
and Lyubarskij et al. No external forcing is applied. 

The equations are first discretized in space by repre- 
senting the variables E and i? on a set of 200 logarithmi- 
cally spaced orbits. The derivatives are represented by sim- 
ple, centred finite differences and the scheme conserves mass 
exactly. The resulting set of temporal ODEs is solved using 
a fifth-order Runge-Kutta method with adaptive step-size. 
The time-step adjusts automatically to ensure accuracy and 
stability of the integration. 

Before the start of the run, the dimensionless coeffi- 
cients (Qi, Q2, Zi/e^'^ , Za/e"^) are evaluated on a rectangu- 
lar grid in the space (e, Ae', Xecj') by solving the dimension- 
less ODEs of Section 4.4 in a once-for-all calculation for the 
chosen parameter values. During the run, the coefficients are 
then evaluated rapidly by trilinear interpolation on the grid. 

Fig. 2 shows the evolution of e and uj over a time inter- 
val loot,/, where t,y = /v is the viscous time-scale at the 
inner boundary, and u = Qf2(X/E) is the vertically averaged 
kinematic viscosity in a circular disc. During this period 
the surface density exhibits negligible evolution. However, 
there is a rapid twisting of the disc because of differential 
precession caused by the radial pressure gradient, an effect 
explained in Section 2.2. After a transient phase, the eccen- 
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^ 0.2 




fundamental set of one-dimensional equations that describe 
how the mass, angular momentum and eccentricity vector of 
a thin disc evolve a s a result of internal stresses and external 
forcing (equations E1S-G20). The analysis is asymptotically 



Figure 2. Top: Evolution of the eccentricity in an initially uni- 
formly eccentric disc. The curves are ordered from top to bottom, 
and are separated by time intervals of lQt,y. Bottom: Evolution 
of the longitude of periastron (radians). The curves are ordered 
from bottom to top. 



tricity settles into a twisted mode that decays exponentially 
in time and precesses slowly in a prograde sense. 

Owing to the monotonic decay of the eccentricity, non- 
linear terms are active only in the earliest stages of the evo- 
lution. However, non- linear effects will be critical in other 
applications, such as determining the outcome of the eccen- 
tric instability in superhump binaries. 



7 DISCUSSION 

In this paper a comprehensive theory of eccentric accre- 
tion discs has been presented. Starting from the basic fluid- 
dynamical equations in three dimensions, I have derived the 



exact in the limit of a thin disc, and allows for slowly varying 
eccentricities of arbitrary magnitude. 

These equations are generally valid and therefore of fun- 
damental interest. They are the equivalent of the Gauss 
perturbation equations for a continuous disc. Previously, 
Lyubarskij et al. (1994) succeeded in deriving a related set of 
equations for an eccentric disc by considering the conserva- 
tion of mass, angular momentum and energy. However, their 
analysis is restricted to the case in which the ellipses are all 
aligned and do not precess. Their method works in this case 
because a knowledge of the angular momentum and energy 
of an orbiting body is sufficient to determine its semi-latus 
rectum and eccentricity, but not its longitude of periastron. 
A closed system of equations is obtained only if the ellipses 
are artificially constrained not to precess. In reality, such 
precession is inevitable and the evolution of the longitude 
of periastron must be determined from a full analysis of the 
horizontal components of the equation of motion. This leads 
to an equation for the eccentricity vector, or complex eccen- 
tricity, which is not in conservative form. 

The second achievement of this paper is the explicit de- 
velopment of the equations in the case of a specific stress 
model which, it is hoped, gives a fair representation of the 
turbulent stress in an accretion disc. To obtain the coeffi- 
cients in the evolutionary equations requires a solution of 
the non-linear PDEs that govern the azimuthal and vertical 
structure of the disc. It also requires an understanding of the 
relation between the turbulent stress tensor and the velocity 
gradient tensor. The simplest plausible relation, adopted in 
almost all theoretical work on accretion discs, is an effective 
viscosity model in which an instantaneous linear relation is 
assumed, and the equation of motion therefore reduces to 
the Navier-Stokes equation. In this paper I have introduced 
a Maxwellian viscoelastic model of the turbulent stress in an 
accretion disc. This generalizes the conventional alpha vis- 
cosity model to account for the non-zero relaxation time of 
the turbulence, and is physically motivated by a considera- 
tion of the nature of MHD turbulence. The PDEs governing 
the azimuthal and vertical structure of the disc, including 
the effects of vertical motion, dissipation of energy and ra- 
diative transport, have been reduced exactly to a set of di- 
mensionless ODEs which can be solved numerically to high 
accuracy to yield the coefficients required for the evolution- 
ary equations. This shows that the technique of non-linear 
separation of variables, applied first to warped discs (Ogilvie 
1999, 2000), is not restricted to purely viscous models but 
can incorporate improved representations of the stress as our 
understanding of magnetorotational turbulence develops. 

It has been confirmed that circular discs are usually vis- 
cously unstable to short- wavelength eccentric perturbations, 
as found by Lyubarskij et al. (1994), if the conventional al- 
pha viscosity model is adopted. It has been noted that the 
instability is essentially the same as the viscous overstability 
of axisymmetric modes discovered by Kato (1978). The in- 
stability can be suppressed by introducing a sufficient effec- 
tive bulk viscosity, although the values required may not be 
realistic. More plausibly, the instability can be suppressed by 
allowing for the non-zero relaxation time of the turbulence. 
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even if the bulk viscosity is zero. It has then been shown 
that an initially uniformly eccentric disc does not retain its 
eccentricity over many viscous time-scales, as had been sug- 
gested by Syer & Clarke (1992, 1993) and Lyubarskij et al. 
(1994). These earlier works neglected the differential pre- 
cession caused by slightly non-Keplerian rotation resulting 
from the radial pressure gradient. This leads to twisting of 
the disc, followed by viscous decay of the eccentricity. 

The theory presented here goes considerably beyond 
previous analytical treatments of eccentric discs. It also pro- 
vides a practical numerical scheme that involves only one- 
dimensional equations and from which the fast orbital time- 
scale has been eliminated. This scheme is to be preferred in 
many circumstances to a direct numerical simulation of the 
fluid-dynamical equations. Almost all direct simulations to 
date attempt to represent the Navier-Stokes equation (with 
or without explicit viscosity) in two dimensions. The present 
analysis shows that a two-dimensional treatment of eccentric 
discs may capture many of the correct qualitative features 
but cannot be trusted in detail. Vertical motion is always 
present in eccentric discs and radiative damping can influ- 
ence the evolution of eccentricity. Furthermore, the Navier- 
Stokes equation does not take into account the relaxation 
time of the turbulence, which can be of great importance in 
this context. 

Nevertheless, it would be valuable to make detailed 
comparisons between the present theory and direct simu- 
lations (preferably three-dimensional). The present theory 
cannot be applied reliably to thick discs, nor to situations 
in which the eccentricity varies rapidly in space (i.e. on 
a length-scale comparable to the thickness of the disc) or 
in time (i.e. on a time-scale comparable to the orbital pe- 
riod). Mean-motion resonances are therefore excluded from 
the analysis, which is secular in the sense of celestial me- 
chanics. Nevertheless, the effect of mean-motion resonances 
could be included in the evolutionary equations by adding 
appropriate localized source terms for angular momentum 
and eccentricity. 

The theory developed in this paper has much in com- 
mon with the theory of warped accretion discs (e.g. Pringle 
1992; Ogilvie 1999, 2000). One distinction, noted above, is 
that the equations for an eccentric disc are not all in conser- 
vative form. Another difference is that the theory of warped 
discs is complicated by a resonance caused by the coinci- 
dence of the orbital and epicyclic frequencies in a Keple- 
rian disc. As a result, the behaviour is qualitatively differ- 
ent depending on the relative magnitudes of a and H/R 
(Papaloizou & Lin 1995). Fortunately, no such complication 



arises iiii Lhe caHU of ticceuU'lc dlHca. A couHiHLouL aaympLoLic 
expansiijii uf the fliiid-d^iiamiLal equaLiuns is possible for 
any value of a, and the fractional error in the asymptotic 
approximation, 0{{H/ R)'^), is very small in most applica- 
tions. 

The evolutionary equations should be useful in many 
applications, including understanding the eccentric planet- 
disc interaction and testing theories of quasi-periodic oscil- 
lations in X-ray binaries. In future work the rate of change of 
the complex eccentricity caused by external forcing should 
be evaluated explicitly in the non-linear case for tidal forc- 
ing by a companion object on a circular or eccentric orbit, 
and also for Einstein precession near a black hole. 
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APPENDIX A: GEOMETRICAL QUANTITIES 

The following dimensionless quantities are required for eval- 
uation of the stress coefficients. These expressions generalize 
those of Lyubarskij et al. (1994) to allow for precession of 
the orbits, and also correct a number of errors in that paper. 
In the following, c and s denote cos6 and sin^, respectively. 
Rate of change of radius with semi-latus rectum: 



7?A = 



1 + (e — Ae')c — \eu' s 
(l + ec)2 ■ 

Inverse metric coefficients: 

(l + ec)^(l + 2ec + e^) 



AA 

9 = 



[1 + {e- Xe')c- Xeuj'sY 
(1 + ec)^es 



, A</) _ 

^ ^ l + (e-Ae')c- Aeon's' 
a'/"* = (1 + ec)^ 

Orbital contribution to the velocity divergence: 

(1 + ec)(Ae's — Xeuj'c — Xe^uj') 
1 + (e — Ae')c — XeLu's 

Variation of the surface density around the orbit: 

^ (l-e^)^/^(l + ec) 
1 + {e- Xe')c- Xeu's' 

Derivatives of the angular velocity: 

93 ^ " -^{i + ec) + 2(Ae'c + Aea;'s)(l + ec), 

g4 ~ — 2es(l + ec). 

)Ut 

(1 + ec) 



Orbital contribution to the shear tensor: 

\3 



AA 

S = 



[1 + (e - Ae')c - Xeuj'sf 
X {(1 + ecfes + Ae'(l + ec + e^s^)s 
-Aetj' [c + e(2 + c^) + e^(4 - c^)c + e^] } 



(Al) 

(A2) 

(A3) 
(A4) 

(A5) 

(A6) 

(A7) 
(A8) 



A^s'^"* = 



(1 + ec)^es 



1 + (e — Ae')c — Xeuj's 
i(l + ec) — Ae'c — Aeon's 



(All) 



Covariant components: 



saa 



(^{-ie(l + ee). + Ae'(l + ^ee) 



'Xeuj' 



c--e(l-3c^) 



,2/2 -.2 ''/-I r,2\|\22/2 "1 

-A e cs — A ee (1 — 2c ) + A e tj csj 



\ -1 

A Sa0 



(1 + ec)2 
3 



X I ec + Ae'c + Xeuj's 1 — Aee' 

4 4 2 



\ -2 

A Stij^ = 



(l + ec)2 



(A12) 

(A13) 
(A14) 



(A9) 



As 



\4, 



(1 + ec) 



[1 + (e - Ae')c - Aetj's]^ 



6-^(1 + 4:0-') e-'c 



+Ae' ' 
+AeLj'(l + 2ec + e^)s} 



-ec - 
4 4 

c- ie(l-4c^) + ie^c 



(AlO) 
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